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ABSTRACT 



Aims. Stellar mass distribution in the Andromeda galaxy (M 31) is estimated using optical and near-infrared imaging data. Combining 
the derived stellar mass model with various kinematical data, properties of the dark matter (DM) halo of the galaxy are constrained. 
Methods. SDSS observations through the ugriz filters and the Spitzer imaging at 3.6 microns are used to sample the spectral energy 
distribution (SED) of the galaxy at each imaging pixel. Intrinsic dust extinction effects are taken into account by using far-infrared 
observations. Synthetic SEDs created with different stellar population synthesis models are fitted to the observed SEDs, providing es- 
timates for the stellar mass surface density at each pixel. The stellar mass distribution of the galaxy is described with a 3-dimensional 
model consisting of a nucleus, a bulge, a disc, a young disc and a halo component, each following the Einasto density distribution 
(relations between different functional forms of the Einasto density distribution are given in Appendix |B](. By comparing the stellar 
mass distribution to the observed rotation curve and kinematics of outer globular clusters and satellite galaxies, the DM halo parame- 
ters are estimated. 

Results. Stellar population synthesis models suggest that M31 is dominated by old (<; 7Gyr) stars throughout the galaxy, with the 
lower limit for the stellar mass-to-light ratios M/L r i 4 Af /L . The upper limit M/L r Sj 6M G /L Q is given by the rotation curve of 
the galaxy. The total stellar mass is (10-15) • 10 10 M o , 30% of which is in the bulge and 56% in the disc. None of the tested DM 
distribution models (Einasto, NFW, Moore, Burkert) can be falsified on the basis of the stellar matter distribution and the rotation 
curve of the galaxy. The virial mass M200 of the DM halo is (0.8-1.1) • 10 12 M G and the virial radius is R20Q =189-213 kpc, depending 
on the DM distribution. For the Einasto profile, the average density of the DM halo within the central 10 pc is 16-61 M Q pc~ 3 (0.6- 
2.3 TeV/e 2 cirT 3 ), depending on the stellar mass model. The central density of the DM halo is comparable to that of nearby dwarf 
galaxies, low-surface-brightness galaxies and distant massive disc galaxies, thus the evolution of central DM halo properties seems to 
be regulated by similar processes for a broad range of halo masses, environments, and cosmological epochs. 

Key words, galaxies: individual: Andromeda (M31) - galaxies structure - galaxies: fundamental parameters - galaxies: kinematics 
and dynamics - galaxies: haloes - cosmology: dark matter 



1. Introduction 

Due to its proximity and size, our nearest large neighbour galaxy 
M 3 1 offers a unique and attractive opportunity to study stellar 
populations and galactic structure in detail. It has been a source 
for groundbreaking discoveries ever since the secure acknowl- 
... edgem ent of this nebula as an extragalactic object by Ernst Opik 
> ■ (119221) . 

Self-consistent treatment of the available photometry and 
5_i 1 kinematical data of M31 enabled the construction of sophis- 
C$ , ticat ed multi-compone nt galactic models already decades ago 
(e.g. lTenies et alj["l994l) . Considering the huge volume and high 
detail of observational information available nowadays, complex 
mass models of M31 offer a promising opportunity for cast- 
ing light on one of the most puzzling problems in astrophysics 
and cosmology: the nature and properties of dark matter (DM) 
haloes. By now, the increasing scope of observations has enabled 
stretching mass models much further than the extent of gas disc 
rotati on curves, providing new clues about DM halo parame- 
ters dGeehan etaljr2 006: Sei gar et al.ll2008i: IChemin et alj [2009: 
ICorbelli et al.ll2010h ~ 



pected flux of high-energy particles resulti ng from decaying or 
annihilating DM in extragalactic sources (Bert one et al.l 12007: 



Hiitsi et al. 1201 Oh . particles from more concentrated extragalac- 
tic objects might be detectable as an enhancement of the Galactic 
signal within certain apertures. By comparing the assumed DM 
distribution in M 3 1 to the data collected with the diverse arsenal 
of ground-based and space-borne detectors of high-energy par- 
ticles, some constraints on the energy spectrum of DM particles 



On the other hand, particle physics instrumentation has 
reached a level at which it can provide some hints about DM. 
Although the diffuse Galactic background likely exceeds the ex- 



have alread^^e^nj^ (ejyAharonian et al. 2003; Lavall e et al 
20061 iBovarskv etal] 120081: iDugger et all 12010b Iwatson et al 
l2012h . 

For more stringent constraints, not only more capable detec- 
tors are needed but also a better understanding of the properties 
of the source DM haloes. As bizarre as it seems, the derivation 
of the detailed mass distribution of the Andromeda galaxy was 
limited by the lack of suitable optical imaging up to recent times. 

Although visible even to the naked eye, its span over 
four degrees on the celestial sphere makes Andromeda a 
real challenge to observe with a usual scientific telescope. 
Thus it is only very recently that observations covering 
the galaxy with deep wide-field CCD imaging have started 
to appe ar: a dedicated scan within the Sloan Digital Sky 
Survey dYork et al.| [2000). the Ca nada-France-Hawaii teles cope 
Megacam programme PAndAS dMcConnachie et afl r2009). the 
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Pan-STARRS telescope project PAndrome da dLee et al.1 20T2I 
Combined with the space-based ultraviolet (Thilkeret al. 2005 
near- (Barmbv et al. 2006) and far-infrared dGordon e t al. 2006; 
iFritz et al.ll201 ll) observations, these data now provide an un- 
precedented panchromatic view of a galaxy at a resolution of a 
few parsecs, allowing for the derivation of detailed properties of 
the stellar populations. 

In this paper we estimate the stellar and DM distribution of 
M31 using the SDSS and Spitzer 3.6-micron imaging to con- 
strain the properties of the stellar populations. We construct a 
mass distribution model of the galaxy in correspondence with 
the latest kinematical data from the literature, giving estimates 
for the DM halo properties and the related uncertainties. 

We have taken the inclination angle of M31 to be 
77.5° (IWalterbos & Kennicutl[l988t |de Vaucouleurs et al. 1991) 
and t he distance to the galaxy 785 kpc (McConn achie et al.l 
120051) . yielding the scale 228 pc/arcmin. Throughout the paper, 
luminosities are presented in AB-magnitud es and are c o rrecte d 
for the Galactic extinction according to Temp el et al.l (1201 ll) . 
where e xtinction correspondi ng to the Sloan filters is derived 
from the lSchlegel et alJ (U998) estimates and the Galactic extinc- 
tion law by linear interpolat ion. The absolute solar luminosity 
for each filter is taken from lBlanton & Roweisl d2007l) . 

2. Observational SEDs 

In an ideal case, one would study stellar populations using all 
the available photometric data to sample the spectral energy dis- 
tribution (SED) of a galaxy throughout the full electromagnetic 
spectrum. However, we have limited ourselves here to the optical 
and near-infrared section of the spectrum, since its interpretation 
with synthetic stellar population models is most straightforward. 
Also the stellar mass is best traced by this wavelength domain. 

For deriving the observed SEDs of M31, we relied on the 
Sloan Digital Sky Survey (SDSS) observations through the ugriz 
filters and the Spitzer Space Telescope IRAC camera imaging 
at 3.6 microns. For our purposes, these observations provide a 
sufficiently wide and deep coverage and the calibration of the 
data is relatively well-understood. 

The basic steps of the SDSS imag e processing and m osaic - 
ing used here have been introduced in Tem pel et alJ (1201 ll) . The 
intrinsic absorption of the galaxy has been t aken into account by 
apply ing the dust disc model developed in Tem pel et al.1 $2010, 
1201 ll) on the basis of the far-infrared flux distribution as mea- 
sured by the Spitzer MIPS camera. We have used the resultant 
absorption-free SDSS images for recovering the total starlight 
along sight-lines affected by the dust disc. The final SDSS mo- 
saic was resampled to 3.96 arcsec px~' . 

To reconstruct the near-infrared view of M 3 1 we retrieved 
the pipeline-calibrated and -processed (post-BCD) Spitzer im- 
ages from the NASA/IPAC Infrared Science Archive. Exposures 
severely suffering from cosmic rays were omitted. A mosaic im- 
age was created using pointing information in the image head- 
ers. The final mosaic was resampled to the same pixel scale as 
the SDSS images. 

The spatial resolution of the dust absorption model was lim- 
ited by the point-spread-function (PSF) of the Spitzer MIPS 
camera at 160 microns, thus we convolved the SDSS and Spitzer 
images with the same PSF. This step improved also the signal- 
to-noise ratio in the outer regions of the galaxy and removed 
negative pixel values resulting from the noise after sky removal. 

The images were matched with each other geometrically to 
sub-pixel accuracy. Foreground stars, satellite galaxies and back- 
ground objects were masked with a common mask for each filter, 



constructed using SExtractor dBertin & Arnoutsf l996) and man- 
ual region placing. Additionally, the noisy edges of the Spitzer 
images were masked. Now the SED could be directly derived by 
calibrating the intensity within each pixel to standard flux units. 

3. Model SEDs 

A wide variety of synthetic stellar population models is avail- 
able for interpreting the observed spectral energy distributions. 
In order to address possible degeneracies and uncertainties of 
such models, we have used three models to reproduce the ob- 
servational SEDs. The chosen models follow significantly dif- 
ferent approaches for generating the properties of the syn- 
thetic stellar popu l ations: (I) the composite model spectra by 
Blant on & Roweisl (120071. hereafter B07) are composed as a 
mixture of iBruzual & Charlotl (120031) instantaneous-burst stellar 
population models of different ages and me tallicities, and mod - 
els of gas emission from MAPPINGS -III dKewlev et al.ll200ll) : 
("ID lMarastonl (|2005 ; i hereafter M05) models lay a special empha- 
sis on the thermally pulsing asymptotic gia nt branch stars; (III) 
the evolutionary synthesis model GALEV ( Kotulla et al. 2009) 
is the sole model in which the chemical evolution of ga s and 
stars is treated self-consistently. In B07, the Chabrier ( 2003) stel- 
lar initial mass func tion (IMF) wa s used; for M05 and GALEV, 
we have chosen the Kroupal (12001) IMF option. 

B07 provides five composite spectra corresponding to an 
extremely young, an old, and three intermediate populations. 
It is shown in B07 that a linear combination of these spectra 
can adequately describe the spectral energy distribution of most 
of the galaxies. However, this aspect alone does not necessar- 
ily prove that the underlying models are meaningful; rather, it 
demonstrates that the spectra are sufficiently diverse. Thus for 
the M05 and GALEV models, we have followed the general 
observational knowledge about the star formation history and 
metallicity of M 3 1 to tame the age-metallicity degeneracy, as- 
suming that much of the galaxy is dominated by old stars of 
nearly solar metalli city, while the star-fo r ming ring is composed 
of younger stars dBellazzini et al.l l2Q03t ISaraiedini & Jablonkal 
2005 : iBrown et al.l 120061: lOlsen et al.l 120061: ISadia et al.ll2Q10h 
IZou et al]|201 ll) . From the available M05 models, we have se- 
lected single-, instantaneous-burst stellar populations. GALEV 
is used in the chemically consistent regime to generate old pop- 
ulations with different star formation histories, with and without 
an additional starburst having occurred billion years ago to 
mimic the star-forming regions. 

For each set of model spectra, we sought linear combinations 
of up to 5 spectra to represent the SED of M 3 1 within each pixel 
according to the formula 



obs 



(i) 



where /(/l) bs is the observed SED, /(/i), are the model SEDs 
per unit mass and m, are the corresponding weights of the model 
SEDs. In this formalism, nij effectively measure the mass of each 
model stellar population within a given pixel. 

The spectra with a non-negligible contribution to the inte- 
gral SED of the galaxy are listed in Table \T\ together with other 
relevant information. An illustration of the SED fitting within 
random pixels from the bulge and disc regions of the galaxy is 
presented in Fig. [1] Here, the observed flux through each of the 
six filters is overlaid with the best-fitting linear combination of 
the B07 composite spectra. 

Following Eq. ([TJ, the SED-fitting process yielded the stellar 
mass of each model population within each imaging pixel. For 
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Table 1. Synthetic stellar populations used for SED fitting. 
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Fig. 1. Examples of the observed (large circles) and modelled 
(lines) SED for a random pixel in the bulge region {upper panel) 
and in the young disc region (lower panel). The sizes of the 
datapoints indicate the photometric uncertainties of each mea- 
surement. The model values corresponding to each filter are also 
shown (small datapoints). In most pixels, the reddest model pop- 
ulation (B07-4) alone provides a good representation of the ob- 
served SED. In the young disc regions, the stellar populations 
are more diverse: in the lower panel, the B07 model populations 
1, 4 and 5 contribute 24.44%, 75.52%, and 0.04% of the mass, 
respectively. The corresponding SEDs are weighted according to 
the mass fraction. In this plot, the observed SEDs and the sum 
of the weighted model spectra are normalised per 1 M Q at a dis- 
tance of 10 pc. 



every population synthesis model, the contribution of different 
synthetic populations to the integral mass of the galaxy is pre- 
sented in the last column of Table [TJ It is seen that for each set 
of model spectra, the reddest spectrum (corresponding to an old 
population with near-solar metallicity) dominates all across the 
galaxy and dictates its population properties. It is only within the 
star-forming ring that other spectra provide detectable contribu- 
tion. 

The mass-to-light ratios and thus the masses of the stellar 
components predicted by different population synthesis mod- 
els are remarkably different, as indicated in Table [TJ M05 and 
GALEV give much more massive stellar populations than B07. 
The scatter of the mass-to-light ratios of different models re- 
sults from differences in the modelling approach and hence re- 
flects the uncertainties of stellar mass estimations in general. 
Although the B07 model spectra provide the most precise de- 
scription of the actual SEDs, we have no grounds to state that 
the B07 model represents the actual populations better than the 
M05 and GALEV models. Thus in the following, we have con- 
sidered this mass scatter as an uncertainty of the final mass esti- 
mates. We used the B07 model for the lowest limit and planned 
initially to use the other models for the upper limit of the stellar 
mass. As shown below, however, the masses suggested by the 
M05 and GALEV models are contradicting the rotation curve 
measurements, which set a more rigid upper limit for the stellar 
masses. Nevertheless, we discourage the reader to use this result 
alone forjudging the M05 and GALEV models in general. We 





A crp 


L rc / n J 


Midi 


Mm 
L, 


Mm 
Li 


Fract 


(i) 


[Gyr] 
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(3) 


r Mo I 
L Is 1 

(4) 


r Mo l 
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(6) 


(7) 


B07-1 


0.7 


0.40 


0.76 


0.78 


0.72 


0.014 


B07-3 


0.4-1 


0.05 


0.47 


0.50 


0.56 


0.003 


B07-4 


7-12 


0.03 


5.05 


3.87 


3.12 


0.983 


M05-1 


1 


0.00 


1.11 


1.00 


0.85 


0.008 


M05-2 


2 


0.00 


2.18 


1.70 


1.43 


0.002 


M05-3 


4 


0.00 


3.99 


3.03 


2.56 


0.214 


M05-4 


12 


0.00 


11.6 


8.08 


6.47 


0.767 


M05-5 


12 


-0.33 


9.00 


6.60 


5.37 


0.009 


GALEV- 1 


1, 10 


0.04 


2.88 


3.14 


2.92 


0.004 


GALEV-2 


2, 11 


0.07 


4.35 


4.13 


3.65 


0.011 


GALEV- 3 


4, 13 


0.09 


7.58 


6.20 


5.23 


0.089 


GALEV-4 


12 


0.12 


4.63 


4.55 


4.05 


0.015 


GALEV-5 


12 


0.18 


10.9 


8.33 


6.86 


0.881 



Notes. The columns contain the following: (1) stellar population 
model; for B07 models the number is as in the original paper; (2) 
approximate age of the dominant star-formation epoch(s); (3) average 
metallicity of the stars; (4)-(6) mass-to-light ratio in the gri filters; (7) 
total stellar mass fraction in M 3 1 of the corresponding stellar popula- 
tion. 
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Fig. 2. Stellar mass-density map of M 3 1 . The ellipses enclose 
50%, 75%, 90%, and 95% of the total mass, respectively. 



can just conclude that with our currently used configuration of 
the stellar population details (i.e. the IMF, metallicity, and star 
formation history), these models overestimate the realistic stel- 
lar mass of M31. 



4. Stellar mass distribution 

In the previous section, the stellar mass corresponding to each 
model spectrum was derived for each imaging pixel. This gives 
the 2-dimensional stellar mass distribution in M 31, as presented 
in Fig. [2] for the B07 model. The distribution appears feature- 
less and regular, resulting from the intensive smoothing with the 
PSF of the Spitzer 160-micron imaging, but it indicates also 
that the galaxy is generally undisturbed and that the intrinsic 
dust-absorption effects have been removed in proper proportions 
along different lines of sight. 

The actual spatial distribution of stellar matter can be split 
into the contributions of different galactic components. We mea- 
sured the elliptically averaged stellar mass distribution from 
Fig.|2]and approximated it as a superposition of the stellar com- 
ponents: a nucleus, a bulge, a disc, and a halo. In addition, the 
ring-like star-forming region was taken as a separate component 
from the disc and is referred here as the young disc. Assuming 



A. Tamm et al.: Visible and dark mass of M 31 



10 u 



5 10 2 



i 10"* 



10" 



3.01 



0.1 



10 



Bulge 



(kpc) 



Disc 

Stellar halo 



Stellar mass profile 
Model profile 
Components 
B07-1 



Young disc: 



10 



100 



1000 



(arcsec) 



Fig. 3. The mass-density distribution of the galaxy, averaged 
along elliptical iso-density contours, as inferred from the B07 
model (thick grey line; its thickness indicates deviations along 
each ellipse), the model profile (solid line) and the contributions 
of the individual stellar components (dashed lines) to the model 
profile. The contribution of the model population B07-1 to the 
mass distribution is also shown; it is closely traced by the young 
disc component of the model. 



each component to be an ellipsoid of rotational symmetry and 
constant axial ratio q, we used the Einasto law 



p(a) = p c exp < -d N 



l/N 



(2) 



to describe the density distribution of a component. Here dis- 
tance from the centre a - yV 2 + z 2 /q 2 , where r and z are the 
two cylindrical coordinates; cIn is a function of N, such that p c 
becomes the density at distance a c , which defines the volume 
containing half of the total mass of the component. The deriva- 
tion of dfi is presented in Appendix [B] Being mathematically 
identical to the Sersic's R 1 ^" model but fitted to the space den- 
sity, Eq. ([2} provides a sufficiently flexible distribution law for 
describing relaxed galactic components. 

For a more accurate description of the star-forming region, 
the spatial density distribution of the young disc is assumed to 
have a toroidal form, approximated as a superposition of a pos- 
itive and a negative density component, both following Eq. (0 
(see the last paragraph of AppendixlBlfor more details). 

The structural parameters of all the stellar components have 
been adopted from a previous analysis dTempel et alJl201 Th . also 
based on the SDSS data. In the referred paper, the Einasto law 
was expressed with respect to the harmonic mean radius in- 
stead of a c . Relations between different functional forms of the 
Einasto distribution are presented in AppendixlBl 

Using the previously derived structural parameters and leav- 
ing the component masses as free variables, we fitted these com- 
ponents to the elliptically averaged stellar mass distribution. The 
lower limit estimates of the masses of the stellar components, de- 
rived from the B07 model, are presented in Table|2]together with 
other related parameters. The upper mass limits are constrained 
by the rotation curve and are 1 .5 times higher for each compo- 
nent (see below). The corresponding mass distribution of each 
stellar component as well as the total stellar mass distribution of 
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Fig. 4. Degeneracies between the masses of the different stellar 
components, shown as likelihood contours of various combina- 
tions of component masses. The final B07 model parameters are 
shown with dots. 

the galaxy are shown in Fig. [3] for the B07 model. To illustrate 
the correspondence between the young disc model component 
and the first model spectrum (B07-1), the contribution of the lat- 
ter to the overall mass distribution is also shown. 

It is natural to suspect that a four-component fit to the galaxy 
stellar mass distribution has to be degenerate to some extent. 
We tested the uniqueness of the model by varying the masses 
of the components and calculating the deviation of each resul- 
tant model fro m the observations using the Bayesian inte rface 
tool MultiNest (iFeroz & Hobsodl2008f iFeroz et al. 2009). The 
results of the degeneracy analysis are presented in Fig. [4] in- 
dicating the likelihood of different combinations of component 
masses. Quite expectedly, the most securely determined com- 
ponent is the bulge and the most unreliable mass estimates are 
for the young disc and halo, both being degenerate with the disc 
mass to some extent. The degeneracies would be higher, if all the 
component parameters were set free. A more conservative two- 
component (bulge + disc) model is described in AppendixlAl 

5. Dynamics and dark matter distribution 

The structural model and masses of the stellar components de- 
rived in Sect, fallowed us to calculate the gravitational potential 
of stellar matter in the galaxy. The gravitational potential of a 
galaxy is also traced by the rotation curve. To match the calcu- 
lated rotation curve with the observed one, the contributions of 
gas and dark matter (DM) have to be added to the stellar mass 
model of the galaxy. 

The contribution of gas to the potential of the galaxy is mod- 
est, thus a precise description of the gas distribution in the model 
is not essential. We have assumed that the distribution of the 
gas disc coincides with that of the young disc simply by rais- 
ing the mass of the young disc b y 6 • 10 9 Mp, which is the ap- 
proxi mate sum of the molecu lar (Dame etal]|1993l:lNieten et alj 
| 2006l) and the neutra l dCarignan et al.ll2006tlChemin et al.ll2 009; 
Cor beili et al.ll2010T) gas mass estimates. Note that the molec- 
ular gas content is rather low in M31, in fact even lower than 
the differences between the neutral gas mass estimates made by 
different authors. 

We have considered various functional forms of DM den- 
sity distribution while incorporating a DM halo component into 
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Table 2. Parameters of stellar components. 



Component 




a c 


<? 


N 


d N 


Pc 


M '' 


M/L g b 


M/L," 


Ml Lj 1 ' 




[kpc] 


[kpc] 








[M Q pc- 3 ] 


[10'°M G ] 


[M /L G ] 


[M /L ] 


[M Q /L ] 


Nucleus 


0.01 


0.0234 


0.99 


4.0 


11.668 


1.713 • 10 u 


0.008 


4.44 


3.20 


2.35 


Bulge 


0.63 


1.155 


0.72 


2.7 


7.769 


9.201 ■ 10~' 


3.1 


5.34 


4.08 


3.01 


Disc 


7.7 


10.67 


0.17 


1.2 


3.273 


1.307 • 10~ 2 


5.6 


5.23 


3.92 


2.92 


Young disc 1 ' 


10.0 


11.83 


0.01 


0.2 


0.316 


1.179 ■ 10~ 2 


0.1 


1.23 


1.12 


0.88 


Stellar halo 


6.3 


12.22 


0.50 


3.0 


8.669 


4.459 ■ 10~ 4 


1.3 


6.19 


4.48 


3.25 



Notes. Structural parameters do, q, N, and galaxy luminosities in SDSS filters (L g , L r , Li) are taken from lTempel et al.l d201 lh . Component masses 
M comp are derived in this paper. (n) Harmonic mean radius (see Appendix |BV {b) Masses and mass-light-ratios corresponding to the B07 model, 
i.e. the lower limits; the upper limits (from the maximum-stellar model) are 1.5 times higher in each case. {c) The structural parameters and p c are 
given for the positive component. In the dynamical models, the gas mass 6 • 1O 9 M is added to the young disc. 



the model galaxy. From observations of the dynamics of galax- 
ies, distributions with a nearly constant inner density (and there- 
fore, "isotherm al" or "cored" haloes) have been derived, e.g. by 
lBurkerd(fl995l) : 



PBurkert(r) = 



( 1 + t)[ 1+ (t) 2 ]' 



(3) 



On the other hand, N-body simulations suggest steeply rising 
DM density toward s the centre (therefore "cuspy" haloes, e.g. 
iMoore etal l dT999h : 



PMoore('") 



Pc 



and lNavarro et all (1 19971) (hereafter NFW): 

Pc 

PNFWO) = — -J. 



(4) 



(5) 



In these equations, p c is a density scale parameter. 

More recently, it has been found that the Einasto distri- 
bution Eq. (f2]i match es with the simulated DM haloes over a 
wider range of rad ii dMerritt et all l2006t iNavarro et all 1201 Ot 
IChemin et al.ll20l"il) . and is gaining popularity for various ap- 
plications. 

Each of the four DM distributions was used in combination 
with the stellar components as determined in Sect.|4]to calculate 
the gravitational potential of the galactic model and the corre- 
sponding rotation curve. 

The model rotation curve was fitted to the observed ro- 
tation curve, composed of two Hi datasets f rom the liter- 
ature: the Westerbork telescope observations (Corbelli et al. 
120101) and data from the Effelsberg and Green Bank telescopes 
(ICarignan et al.1 120061) . We did not attempt to include observa- 
tions from the inner parts of the galaxy, where the dynamics of 
gas clouds is too much affected by non-circular motions, leading 
to difficulties in interpreting the data. 

In addition to the gas rotation curves, we used circular veloc- 
ities calculated from the measurements of the motions of glob- 
ular clusters, satellite galaxies, and stellar streams (Table 0), al- 
lowing us to trace the gravitational potential of M 3 1 out to a 
projected distance of more than 500 kpc from the centre. 

For fitting the model rotation curve to the observed one, the 
DM halo parameters were left free, while the masses of the stel- 
lar components were kept fixed. During the first runs of the fit- 
ting, the Einasto shape parameter N was allowed to vary freely, 
which lead to a wide variation of its value. To reduce degrees 
of freedom, we app lied a fixed value at N = 6.0 according to 
iMerritt et alJ {2006) and Nav arro et al.l d2010h . 
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Fig. 5. Upper panel: the observed rotation curve (data points 
with error bars) overplotted with the model (solid line). 
Contributions of each component are also shown (dashed lines). 
The model corresponds to the B07 stellar mass estimates and the 
Einasto distribution for the DM density. Lower panel: the same 
stellar model with four different DM density distributions. For 
clarity, only the total rotation curves and the DM contributions 
are shown. 



As shown in Sect. [3] the stellar masses yielded by different 
stellar population synthesis models vary by a factor of two. We 
first considered two cases: the lowest-mass model with the B07 
mass estimates and the highest-mass model with the other mass 
estimates. In the latter case, however, the stellar mass becomes 
too high, raising the model rotation curve above the observed 
values at distances 10-20 kpc from the centre even without the 
inclusion of a dark matter component (a DM halo would still 
be needed though to gain a match with the outer rotation curve 
observations). Thus we had to abandon the idea of determining 
the upper limits of the masses of the stellar components by using 
stellar population synthesis models. 

In cases when stellar masses cannot be estimated indepen- 
dently, the maximum-disc approach is often followed, i.e. first 
a pure disc is fitted to the observed rotation curve with as high 
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Table 4. DM halo parameters for various distribution functions. 
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1.40 ±0.27 
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44 
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1.54 • 1Q- 2 


37.95 



Notes. (n) Parameter N has been taken 6.0, yielding d N = 17.668. Spherical symmetry is assumed by taking q ■■ 
(A) Dark matter parameters for the maximum-stellar model. 



1 in Eq. l[2j, in which case r c = a c . 



Table 3. Enclosed mass estimates and the corresponding circular 
velocities at large galactocentric radii of M 3 1 . 
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mass as possible and the other components are added thereafter 
in the required proportions. In the case of M 31, the disc mass- 
to-light ratio is very close to the bulge one because of their simi- 
lar ages and metallicities. Therefore, instead of a maximum-disc 
approach, we applied a maximum-stellar model, conserving the 
relative values of the mass-to-light ratios of the stellar compo- 
nents as determined with stellar population synthesis models, but 
multiplying them with a common constant. Without the inclu- 
sion of a DM halo, the rotation curve allowed the multiplication 
of the B07 stellar masses bye 1 .5 at maximum. Henceforth, we 
are using the corresponding model (together with a minimally 
required DM halo) as an upper limit of the stellar masses and 
refer to it as the maximum stellar model. 

Parameters of the best-fitting DM models, corresponding to 
the different distribution functions and the B07 stellar mass es- 
timates are presented in Table |4] For the Einasto DM distribu- 
tion, parameters corresponding to the maximum-stellar model 
are also given. 

The upper panel of Fig. [5] presents the observed rotation 
curve, over-plotted with the curve derived from the B07 stel- 
lar masses and the Einasto DM halo. Contributions of each stel- 
lar component and the DM halo are also shown. In the lower 
panel, model rotation curves for all the four DM models are pre- 
sented. It is seen that within the observed range of the rotation 
curve, differences between different DM profiles are negligible. 
From 7 kpc inwards along the major axis, outside the range of 




flproj [kpc] 

Fig. 6. Outer rotation curve observations and models (upper 
panel), calculated for the B07 stellar masses, and the corre- 
sponding cumulative mass {lower panel). With the exception of 
the Burkert distribution, all DM models fit well to the observa- 
tions. 



observations, the model with the Burkert DM profile (and to a 
lesser extent, also the model with the NFW DM profile) starts 
to deviate from the other models. In Fig. [6] the outer rotation 
curve (upper panel) and the corresponding enclosed mass (lower 
panel) are shown together with the model curves. Again, all the 
DM distribution models match the observations within the er- 
rorbars, except for the Burkert distribution, which produces a 
slightly lighter DM halo, missing a few outer datapoints. 

The model rotation curves for the maximum-stellar model 
are shown in Fig. [7] Now the fit is somewhat worse than for 
the B07 model, especially at the innermost observational data- 
points, securing that the maximum-stellar model indeed provides 
the very upper limits for the stellar masses. 

As shown above, it is not possible to prefer any of the 
given DM distribution models on the basis of the data on M 3 1 . 
Furthermore, in each case, the derived characteristic radii and 
densities of the DM haloes are very degenerate, as indicated in 
Fig- IS] a significant increase of the characteristic radius can eas- 
ily be compensated by lowering the characteristic density and 
vice versa. In this plot, the virial mass M200, defined as the 
mass within a sphere of mean density 200 times the cosmologi- 
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Fig. 7. The observed rotation curve together with the maximum- 
stellar model, in which the stellar masses are 1.5 times higher 
than in the B07 model. 



cal critical densitjQ, (e.g. iNavarro et al] 1201 Oh . is also shown in 
colour coding for each DM model. Interestingly, despite the un- 
certainty of the DM density distribution, the virial mass is well 
constrained, regardless of the chosen DM distribution model. 
The same statement holds in Fig.|9j where the Einasto DM halo 
parameter likelihoods for the two stellar mass models are com- 
pared. The virial mass is actually quite firmly established by 
the outer "test particles" of the dynamics and is almost inde- 
pendent of the stellar model of the galaxy. For the "cuspy" DM 
profiles (Einasto, NFW, Moore), the derived virial mass is ( 1 .04- 
1.13) • 1O 12 M and the corresponding virial radius 207-213 kpc. 
For the "cored" Burkert profile, the values are 0.8 ■ 1O 12 M and 
189 kpc, respectively. In the case of the Einasto DM distribu- 
tion, a sphere extending to 10 pc from the centre contains 16- 
61 M pc~ 3 (0.6-2.3 TeV/c 2 cm 4 ) of DM on an average. 

6. Discussion and comparison to previous models 

Let us compare our mass models of M31 to some other re- 
cently constructed models. In Table[5] mass estimates suggested 
by our models for t he bulge, disc, and DM halo are compare d 
to the estima t es by iGeehan et alJ d2006h. fSeiear et alJ ([2008), 
IChemin et ail d2009l) . and lCorbelli et all d2010t) . For a better un- 
derstanding of the compatibility of these models, we will briefly 
summarise the principal properties and differences of these mod- 
els below. 

The model derived bv lGeehan et alJ d2006l) consists of a cen- 
tral supermassive black hole, a bulge, a disc, and a DM halo. 
Its stellar components are determined using various luminosity 
measurements in V, R, and r filters out to the (projected) dis- 
tance of 25 kpc along the major axis, and the kinematics is based 
on a composite rotation curve. Mass-to-light ratios of the stel- 
lar components are treated as free parameters. The total mass is 
constrained by data on the motions of outer plan etary nebulae, 
globular clusters and satellite galaxies. In Table |5j Geehan et al. 
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r c (kpc) 
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Virial Mass (10 10 M sun ) 

Fig. 8. Parameter likelihoods for different DM density distribu- 
tions in the case of B07 stellar masses. The virial mass corre- 
sponding to each parameter combination is shown in colour cod- 
ing; 90-10 10 M o , 110-10 10 M o , and 130 -10 10 M o levels are indi- 
cated with solid contours. For the Einasto, NFW, and Moore DM 
profile, the virial mass is almost the same. 



1 Here the critical density is calculated for the Hubble constant value 
o = 71kms~ 1 Mpc~ 1 . 



(2006) values for the best-fit model and the maximum-disc 
model (in brackets) are presented. 

Seiga r et al.l d2008l) constructed a similar black 
hole + bulge + disc + DM halo model. They used the Spitzer 
3.6-micron imaging data and B-R colour profile to determine 
the mass-to-light ratios ; dynamical mass estimators were the 
same or similar as in Geeh an et a I] d2006h . In addition to the 
usual DM profiles, Sei gar et alJ (l2008h considered the case of 
a dark halo that has undergone an adiabatic contraction due to 
the gravitational attraction of the baryonic material. In Table [5] 
the model with adiabatic contraction is presented in brackets; it 
should not be compared directly to the other models. 

IChemin et al. d2009h constructed a black 
hole + bulge + disc + gas model, using their newer H i data 
for constraining the kinematics. For a more accurate description 
of the disc potential, the disc density distribution was settled 
as the residual of the surface brightness distribution after the 
subtraction of the bulge contribution. The small contribution of 
gas mass was considered on the basis of H2 and Hi surveys. In 
Table|5] the "hybrid" model (with the bulge mass is determined 
from stellar velocity dispersions and t he disc mass f r om si mple 
stellar population models) values of IChemin et alJ d2009l) are 
pre sented. 

ICorbelli et a l. (2010) used a bulge + disc + gas model to- 
gether with an optional Burkert/NFW DM halo to fit their H 1 
kinematics data and outer dynamics estimators. The best-fit 
model with equal bulge and disc mass-to-light ratios and the 
NFW dark halo is shown in Ta ble [5] We have rescaled the virial 
mass Mag given bv lCorbelli et all (1201 Oh to M 2 oo- 
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Table 5. Comparison of bulge, disc, and DM halo mass estimates. All masses are in 10 ll) M ( 



Model 






M 200 


Geehan et al. (2006), best-fit (maximum-disc) model 


3.3 


8.4(13.7) 


68 (94) 


Seigar et al. (2008), model without (with) adiabatic contraction 


3.5 (3.5) 


5.8(7.3) 


73 (89) 


Chemin et al. (2009). "hybrid" model 


2.32 


7.1 


100 


Corbelli et al. (2010). NFW model with (M/Lk,,„, = (M/Lh-^ 


3.8 


8.8 


85' 


This work, B07 model 


4.4" 


5.7'' 


113 


This work, maximum-stellar model 


6.6° 


8.6'' 


127 



Notes. (fl) Sum of the bulge and stellar halo masses. (6) Sum of the disc and young disc masses. (c) Recalculated from M 98 in the original paper. 
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Fig. 9. Comparison of DM halo parameters for B07 and max- 
imum stellar mass models. In the upper panel, the virial mass 
corresponding to each parameter combination is shown in colour 
coding; 90- 10 10 M o , 1 10-10 10 M o , and 130 -1O 1O M levels are in- 
dicated with solid contours. In the lower panel, virial radii are 
shown in colour coding, with solid contours tracing 200, 210, 
and 220 kpc. 



In contrast to these four works, our stellar model is based 
on fully 2 -dimensional dust-corrected imaging through 6 filters 
and some more recent dynamical mass estimators. Instead of the 
central black hole, our model includes the nucleus of the galaxy, 
which contributes significantly more to the total mass and the 
model rotation curve of the galaxy. This contribution is neverthe- 
less tiny and has a negligible effect on other model parameters, 
as well as our usag e of an oblate bulge (with an axial ratio 0.8; 
Temp el et al.ll201 II) instead of a spherical one. In Table [5] our 
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LSB-NFW: Coccato el al. (2008) 
LSB-iso: Coccato et al. (2008) 
LSB-iso: de Blok et al. (2008) 
LSB-iso: Kuzio de Naray et al. (2008) 
Dw-NFW: Breddels et al. (2012) 
Dw-NFW: Adams et al. (2012) 
Dw-NFW: Gilmore et al. (2007) 
Dw-iso: Gilmore et al. (2007) 
Dw-iso:Oh etal. (2011) 
Dw-iso: Gentile et al. (2007b) 
Dw-iso: Gentile et al. (2007a) 
Einasto (max-stellar model) 
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Fig. 10. Average DM density inside a given radius, correspond- 
ing to different DM distributions in the case of the B07 stellar 
masses. For the Einasto DM distribution, also the maximum- 
stellar-mass case is plotted. For comparison, central densities of 
some nearby dwarf galaxies and low-surface-brightness galaxies 
are shown. The triangular/quadrangular datapoints are calculated 
assuming the Burkert DM, the circular datapoints correspond to 
the NFW DM. 



B07 model and the maximum-stellar model results are given. 
The actual values probably lie between the estimates of these 
two models. 

Table |5] shows that the bulge mass suggested by our models 
is somewhat higher than in previous models, probably resulting 
from different bulge parameters but also because we have used 
a larger set of observational data to constrain the stellar popu- 
lations. Nevertheless, the bulge dominates the total gravitational 
potential only up to the radius 6-8 kpc and bulge properties have 
little effect on the DM halo parameters. 

As can be seen from Table [5] the DM halo virial mass 
estimates have previously remained mostly below 1-1O 12 M , 
whereas our models suggest slightly higher values, (1.0- 
1.3)-10 12 M G . Once again, the most likely source of differences 
is our usage of a larger collection of observational data on the 
outer dynamics. As shown in Figs.[8]and|9] the virial mass is al- 
most independent of the DM density profile and the stellar mass 
model, being uniquely determined by the outer dynamics of the 
galaxy. 

In Fig. [10] the distribution of the average DM halo density 
within a given radius is shown for each DM density distribu- 
tion model. The added datapoints provide an illustrative com- 
parison with DM haloes of other galaxies, for which the aver- 
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age central density has been r ecently measured more or less reli- 
ably: nearby dwarf galaxies (Gilmore et al. 2007; Gentil e et alJ 
2007aE lOh et al.1 1201 It lAdams et alJ 120121: iBreddels et al l 
20121) and low-surface-brightness galaxies flCoccato et al. 2008; 
de Blok et al1l2008t iKuzio de Narav et alj|2008l) . Some of these 
values are calculated for the modified-isothermal DM distribu- 
tion, others for the NFW distribution. In several cases, both ver- 
sions are presented. These datapoints should be compared to the 
Burkert and NFW profiles of M31, respectively. It is seen that 
the central density of DM haloes varies by a couple of orders of 
magnitude and despite its higher total mass, the DM halo of M 3 1 
cannot be distinguished from an average dwarf or low-surface- 
brightness galaxy in this aspect. Interestingly, the estimate of the 
central density (0.012-0.028)-M H pc ~ 3 of DM haloes of massive 
disc galaxies near redshift z 0.9 dTamm & Teniesl l2005) also 
falls within this range, hinting that the DM halo concentration 
process seems to be restricted rather uniformly over a very wide 
variety of halo masses, environments, and cosmological epochs. 

To conclude our work, it is interesting and also disappointing 
to note that the usage of additional observational data does not 
reduce significantly the uncertainties and scatter of the parame- 
ters of M31 mass distribution models. Our vague understand- 
ing of the evolution of the physical properties of stellar pop- 
ulations restrains the gain from all the gathered observational 
information on the chemical content and formation history of 
a galaxy. Despite the improved imaging and kinematics data, 
we are still unable to confirm or rule out the maximum-disc or 
maximum-baryonic approach in splitting the contributions of lu- 
minous and dark matter to the overall mass distribution. It can 
only be concluded that the bulge mass of M 3 1 probably lies 
within the range (4.4-6.6) -1O IO M and the disc mass within the 
range (5.7-8.6) TO IO M . Nevertheless, M31 provides an excep- 
tional opportunity to estimate the virial mass and the outer dis- 
tribution of a DM halo thanks to the possibility of tracing the 
gravitational potential with various test bodies at large distances 
from the galactic centre. 
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Fig.A.2. The observed inner (upper panel) and outer (lower 
panel) rotation curve compared to the model with two stellar 
components. 



Fig. A.l. Mass profile for two stellar components. For major (up- 
per panel) and minor (lower panel) axis. 

Table A.l. Parameters for the two-component model. 
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the form as originally defined and used by Jaan Einasto ( 1969) 
and thereafter, in which the density distribution is given in re- 
spect to the harmonic mean radius ao- Galaxy component is ap- 
proximated by an ellipsoid of revolution with a constant axial 
ratio q; its spatial density distribution follows the law 



p(a) = p exp 



a 

kao 



l/N 



(B.l) 



where po = hM I (Anqafy is the central density and M is the com- 
ponent mass; a = yV 2 + z 2 /q 2 , where r and z are two cylindrical 
coordinates. The coefficients h and k are normalising parameters, 
dependent on the structure parameter TV. 

A detailed derivation of the constant / ; and k depending on TV 
are given in Appendix B of Tenie s et al.l (119941) . In short, to find 
h and k, the following equality must be satisfied 



Ax. (B.2) 



Appendix A: Simple bulge + disc model of M 31 

For several applications (e.g. for comparing to models of more 
distant galaxies), a model of M 3 1 with two stellar components 
is accurate enough. In Fig. lA.ll we present such a model, con- 
servatively taking N — 4 for the bulge and N = 1 for the disc. 
The other parameters are determined by fitting the model com- 
ponents to the B07 stellar mass density distribution derived in 
Sect. |U they are presented in Table IA.1I The corresponding ro- 
tation curve, including the contribution of the Einasto DM halo 
as derived above, is compared to the observations in Fig. IA.2l 



Appendix B: Relations between different forms of 
Einasto distribution 



In our prev ious papers dTempel & Teniesl l2006t iTempel et al.l 
1201(1120111) . we have used the Einasto distribution function in 



These integrals can be solved analytically, giving 

, r(2/v) , r 2 (3A0 



F(3A0' 



AT 3 (2A0' 



(B.3) 



where F is the (complete) gamma function. 

The advantages of this form are the usage of the harmonic 
mean radius as a good characteristic of the real extent of the 
component, rather independent of the structural parameter N, 
and the directly calculable integral mass (or total luminosity). 

More recent ly, instead of the u sual Navarro-Frenk-White 
(NFW) formula jkavarro et alJI 1997b. the Einasto law has been 
used for DM haloes dNavarro et al.ll2004t iMerritt et al.ll2006l) . In 
these works, the equation is used in the form 



p(a) =p c expl-d, 



1/iV 



(B.4) 
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Fig. B.l. Upper panel shows the mass fraction for a radius a c as 
a function of jV, where dx is found using the Eq. ( IB. 6b . Note the 
high amplification level of the y-axis scale. Lower panel shows 
the relation between harmonic mean radius arj> half-mass radius 
a c , and radius a_2i ^jv is calculated using the Eq. (IB. 5b . 

where djy is a function of N, such that a c defines a volume con- 
taining half of the total mass and p c becomes the density at radius 
a c \ a is the same as in Eq. (IB. lb . This Eq. ( IB.4l i is in fact the same 
as used by J. Einasto in his first paper dEinastoll965l) to describe 
the spatial density distribution of galactic components. 

An integral of Eq. ( IB .4! ) over some volume gives the enclosed 
mass, which is finite and in the spherical case (q — 1), equal to 



M(a) 



a 

- 4 *f 



p(x)x 2 dx = AnNa]p c —y 
"n 



3N, d N | — 

a c 



, (B.5) 



where y is the lower incomplete gamma function. 

Replacing y[3N,d N (a/a c ) 1/N ] with T(3A0 in Eq. IB31 gives 
the total mass of the component. 

The val ue of d\ can be fo und exactly, by solving F(3A0 = 
2y(3N, Q?w). lMerritt et all d2006l) showed that the value of d N can 
be approximated by the expression 



3N- 1/3 + 0.0079/N for N > 0.5. 



(B.6) 



A slightly more exact approx imation for d^ is derived in 
Retana-Montenegro et al. (2012), applicable also for Af t 0.5. 

In the upper panel of Fig. IB. II we show the mass fraction 
within radius a c as a function of N, where dn is found from 
Eq. (IB. 6b . It is seen that for N ^ 0.5, Eq. (IB. 6b gives a good 
approximation. For lower N values, the exact solution should be 
preferred for d^. Throughout this work, we have used and pre- 
sented the exact dn values. 

In addition to Eqs. ( IB. lb and (IB ,4b . the Einasto distribution 
for DM haloes has been used in the form 



p(a) =p_ 2 exp-i — 



a 

a-2 



- 1 



(B.7) 



where a — l/N and a_2 marks the radius where the log- 
arithmic slope of the profile^ equals the isothermal value 



2 In fact, in the original paper Einasto ( 1965) introduced just this kind 
of formula. Instead of p_ 2 he used a logarithmic slope in the solar neigh- 
bourhood. 



-dlnp/dlna = 2. This form of the Einasto law is gaining pop- 
ularity for de scribing DM haloes in A^-body simulations (e.g. 
iNavarro et alj |2010[). as well as applications in astroparticle 
physics (e.g. lHutsi et al]l2010t Qfempel et a p2012l) . 

All forms of the Einasto law (Eqs. IB.U IB. 41 and IB. 7b are 
identical, and strict transductions exist between the parameters: 



a c — aok (dN) N , 

p c = p exp(-d N ), 

a- 2 = a k(2N) N = a c (2N/d N f , 

P-2 = Po exp (-2A0 = p c exp (d N - 2N) . 



(B.8) 
(B.9) 
(B.10) 
(B.ll) 



The lower panel in Fig. IB. II shows the relation between the 
harmonic mean radius a®, half-mass radius a c , and radius a_2- 
It is seen that with increasing Af the profile becomes more ex- 
tended and also the half-mass radius increases, while a^ de- 
creases compared to the harmonic mean radius ao- 

Observations have demonstrated that the stellar disc can have 
a toroidal form in some galaxies, i.e. it does not continue to the 
centre (one example is currently studied galaxy, M 31). One sim- 
ple way to represent a stellar population with the central density 
minimum in the present ed density distribution framework was 
introduced by J. Einasto dEinastoll969l:lEinasto et al.lll980h . The 
spatial density of a disc with a central hole can be expressed as 
a sum of two spheroidal mass distributions 



Pdisc(a) =p+(a)+p_(a), 



(B.12) 



both of which can be approximated with the exponential law 
(IB. lb . Characterising the extent of the minimum with a param- 
eter k = flo_/ao+ and adopting a model disc with a zero density 
at r — and a non-negative density at pdisc(a) > 0, the following 
relations must hold: M_ = -k 2 M+, q_ = q+/K, where k < 1 is 
a parameter that determines the relative size of the hole in the 
centre of the disc. The structural parameters A^_ and A^+ should 
be identical. 
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